/*
** LuaFileSystem
** Copyright Kepler Project 2003 - 2017 (http://keplerproject.github.io/luafilesystem)
*/
//#define EIGEN_USE_MKL_ALL

#include <Eigen/Sparse>
#include <Eigen/SparseCholesky>
#include <Eigen/SparseLU>
#include <log/log.h>

#include <unsupported/Eigen/ArpackSupport>

template<class T>
int solve(Eigen::SparseMatrix<T>& A,const Eigen::Matrix<T, Eigen::Dynamic, 1>& b,Eigen::Matrix<T, Eigen::Dynamic, 1>& x)
{
	typedef typename Eigen::Matrix<T, Eigen::Dynamic, 1> Vector;
	typedef typename Eigen::SparseMatrix<T> Matrix;

	Eigen::SimplicialLDLT<Matrix> solver;
	solver.compute(A);
	if(solver.info()!=Eigen::Success) {
	  // decomposition failed
	  return -1;
	}

	x = solver.solve(b);
	if(solver.info()!=Eigen::Success) {
	  // solving failed
	  return -1;
	}
	return 0;
}

template<class T>
int solve_sp_lu(Eigen::SparseMatrix<T>& A,const Eigen::Matrix<T, Eigen::Dynamic, 1>& b,Eigen::Matrix<T, Eigen::Dynamic, 1>& x)
{
	typedef typename Eigen::Matrix<T, Eigen::Dynamic, 1> Vector;
	typedef typename Eigen::SparseMatrix<T> Matrix;

	Eigen::SparseLU<Matrix> solver;
	solver.compute(A);
	if(solver.info()!=Eigen::Success) {
	  // decomposition failed
	  throw("decomposition failed");
	  return -1;
	}

	x = solver.solve(b);
	if(solver.info()!=Eigen::Success) {
	  // solving failed
	  return -1;
	}
	return 0;
}

template<class T>
int solve_sp_lu(Eigen::SparseMatrix<T>& A,
			 const Eigen::SparseMatrix<T>& B,
			 Eigen::SparseMatrix<T>& X)
{
	typedef typename Eigen::Matrix<T, Eigen::Dynamic, 1> Vector;
	typedef typename Eigen::SparseMatrix<T> Matrix;

	Eigen::SparseLU<Matrix> solver;
	solver.compute(A);
	if (solver.info() != Eigen::Success) {
		if (solver.info() == Eigen::NumericalIssue){
			logs::warning_log(
					"The provided data did not satisfy the prerequisites.");
		}
		if (solver.info() == Eigen::NoConvergence){
			logs::warning_log("Iterative procedure did not converge.");
		}
		if (solver.info() == Eigen::InvalidInput){
			logs::warning_log(
					"The inputs are invalid, or the algorithm has been improperly called.");
		}
		return -1;
	}

	X = solver.solve(B);
	if(solver.info()!=Eigen::Success) {
	  // solving failed
	  return -1;
	}
	return 0;
}





template<class T>
int solve_sp_ldlt(Eigen::SparseMatrix<T>& A,const Eigen::Matrix<T, Eigen::Dynamic, 1>& b,Eigen::Matrix<T, Eigen::Dynamic, 1>& x)
{
	typedef typename Eigen::Matrix<T, Eigen::Dynamic, 1> Vector;
	typedef typename Eigen::SparseMatrix<T> Matrix;

	Eigen::SimplicialLDLT<Matrix,Eigen::Upper> solver;
	solver.compute(A);
	if(solver.info()!=Eigen::Success) {
	  // decomposition failed
	  return -1;
	}

	x = solver.solve(b);
	if(solver.info()!=Eigen::Success) {
	  // solving failed
	  return -1;
	}
	return 0;
}


template<class T>
int solve_sp_cg(Eigen::SparseMatrix<T>& A,const Eigen::Matrix<T, Eigen::Dynamic, 1>& b,Eigen::Matrix<T, Eigen::Dynamic, 1>& x)
{
	typedef typename Eigen::Matrix<T, Eigen::Dynamic, 1> Vector;
	typedef typename Eigen::SparseMatrix<T> Matrix;

	Eigen::ConjugateGradient<Matrix,Eigen::Upper> solver;
	solver.compute(A);
	if(solver.info()!=Eigen::Success) {
	  // decomposition failed
	  return -1;
	}

	x = solver.solve(b);
	if(solver.info()!=Eigen::Success) {
	  // solving failed
	  return -1;
	}
	logs::info_log(boost::format("CG iterations: %1%, error: %2%") % solver.iterations() % solver.error());
	return 0;
}

template<class T,class T2>
int solve_sp_eig(Eigen::SparseMatrix<T>& A,
			  Eigen::Matrix<T2, Eigen::Dynamic, 1>& v,int nbrEigenvalues)
{
	typedef Eigen::SparseMatrix<T> SparseMat;
	typedef Eigen::SimplicialLDLT<SparseMat> SparseChol;
	typedef Eigen::ArpackGeneralizedSelfAdjointEigenSolver <SparseMat, SparseChol> Arpack;
	Arpack solver;
	solver.compute(A, nbrEigenvalues, "SM");

	if(solver.info()!=Eigen::Success) {
	  // solving failed
	  return -1;
	}
	v=solver.eigenvalues();
	//arpack.eigenvectors();
	return 0;
}

template<class T,class T2>
int solve_sp_geig(Eigen::SparseMatrix<T>& A,
				Eigen::SparseMatrix<T>& B,
				Eigen::Matrix<T2, Eigen::Dynamic, 1>& lambda,
				Eigen::Matrix<T2, Eigen::Dynamic,  Eigen::Dynamic>& V,
				int nbrEigenvalues)
{
	typedef Eigen::SparseMatrix<T> SparseMat;
	typedef Eigen::SimplicialLDLT<SparseMat> SparseChol;
	typedef Eigen::ArpackGeneralizedSelfAdjointEigenSolver <SparseMat, SparseChol> Arpack;

	Arpack solver;
	solver.compute(A, B, nbrEigenvalues, "SM");

	if (solver.info() != Eigen::Success) {
		if (solver.info() == Eigen::NumericalIssue){
			logs::warning_log(
					"The provided data did not satisfy the prerequisites.");
		}
		if (solver.info() == Eigen::NoConvergence){
			logs::warning_log("Iterative procedure did not converge.");
		}
		if (solver.info() == Eigen::InvalidInput){
			logs::warning_log(
					"The inputs are invalid, or the algorithm has been improperly called.");
		}
		return -1;
	}
	lambda=solver.eigenvalues();
	V = solver.eigenvectors();
	return 0;
}
